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Quantum chromodynamics (QCD) is the theory of the strong interaction, explaining (for 
] example) the binding of three almost massless quarks into a much heavier proton or neutron 
. and thus most of the mass of the visible Universe. The standard model of particle physics 

_ predicts a QCD-related transition that is relevant for the evolution of the early Universe. At 
i low temperatures, the dominant degrees of freedom are colourless bound states of hadrons (such 
J> [ as protons and pions). However, QCD is asymptotically free, meaning that at high energies or 
O ■ temperatures the interaction gets weaker and weaker [1; 2], causing hadrons to break up. This 
behaviour underlies the predicted cosmological transition between the low-temperature hadronic 
phase and a high-temperature quark gluon plasma phase (for simplicity, we use the word 'phase' 
to characterize regions with different dominant degrees of freedom). Despite enormous theoretical 
effort, the nature of this finite-temperature QCD transition (that is, first-order, second-order or 
analytic crossover) remains ambiguous. Here we determine the nature of the QCD transition 
using computationally demanding lattice calculations for physical quark masses. Susceptibilities 
are extrapolated to vanishing lattice spacing for three physical volumes, the smallest and largest 
of which differ by a factor of five. This ensures that a true transition should result in a dramatic 
increase of the susceptibilities. No such behaviour is observed: our finite-size scaling analysis 
shows that the finite-temperature QCD transition in the hot early Universe was not a real phase 
' transition, but an analytic crossover (involving a rapid change, as opposed to a jump, as the 
temperature varied). As such, it will be difficult to find experimental evidence of this transition 
from astronomical observations. 

During the evolution of the Universe there were particle-physics-related transitions. Although there are 
^ i' strong indications of an inflationary period, we know little about how it affected possible transitions of our 
D . known physical model. To understand the consequences, we need a clear picture about these cosmologically 
relevant transitions. The Standard Model (SM) of particle physics predicts two such transitions. 

One of the SM based transitions occurs at temperatures (T) of a few hundred GeV. This transition is 
responsible for the spontaneous breaking of the electroweak symmetry, which gives the masses of the elementary 
particles. This transition is also related to the electroweak baryon- number violating processes, which had a major 
' influence on the observed baryon-asymmetry of the Universe. Lattice results have shown that the electroweak 
transition in the SM is an analytic cross-over [3; 0; S [sl . 

The second transition occurs at T «200 MeV. It is related to the spontaneous breaking of the chiral symmetry 
of QCD. The nature of the QCD transition affects our understanding of the Universe's evolution (see Ref. [7| 
for example). In a strong first-order phase transition the quark-gluon plasma supercools before bubbles of 
hadron gas are formed. These bubbles grow, collide and merge, during which gravitational waves could be 
produced [8]. Baryon-enriched nuggets could remain between the bubbles, contributing to dark matter. The 
hadronic phase is the initial condition for nucleosynthesis, so inhomogeneities in this phase could have a strong 
effect on nucleosynthesis As the first-order phase transition weakens, these effects become less pronounced. 
Our calculations provide strong evidence that the QCD transition is a crossover and thus the above scenarios 
— and many others — are ruled out. 

We emphasize that extensive experimental work is currently being done with heavy ion collisions to study 
the QCD transition (most recently at the Relativistic Heavy Ion Collider, RHIC). Both for the cosmological 
transition and for RHIC, the net baryon densities are quite small, and so the baryonic chemical potentials {n) 
are much less than the typical hadron masses («45 MeV at RHIC and negligible in the early Universe). A 
calculation at ii=0 is directly applicable for the cosmological transition and most probably also determines the 
nature of the transition at RHIC. Thus we carry out our analysis at /i=0. 

QCD is a generalised version of quantum electrodynamics (QED). The Euclidean Lagrangian with gauge 
coupling g and with a quark mass of m can be written as £=— l/(2(7^)TrF^j^F^i,-|-'!/j7^(9p-|-A^-|-m)'i/;, where 
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Ffj.iy—d^A,y-d,^Af^+[Af^,Ai,]. In electrodynamics the gauge field A^ is a simple real field, whereas in QCD it is 
a 3x3 matrix. Consequently the commutator in i^^jy vanishes for QED, but it does not vanish in QCD. The 
ip fields also have an additional "colour" index in QCD, which runs from 1 to 3. Different types of quarks are 
represented by fermionic fields with different masses. The action S is defined as the four-volume integral of C. 
The basic quantity we determine is the the partition function Z, which is the sum of the Boltzmann factors 
exp(— 5') for all field configurations. Partial derivatives of Z with respect to the masses give rise to the order 
parameters we studied here. 

There are some QCD results and model calculations to determine the order of the transition at fi—O and 
/i^^O for different fermionic contents (compare refs 13; U; 13; 13; 14; l^; TEi, U; 18; l3|)- Unfortunately, none 
of these approaches can give an unambiguous answer for the order of the transition for physical values of the 
quark masses. The only known systematic technique which could give a final answer is lattice QCD. 

Lattice QCD discretises the above Lagrangian on a four-dimensional lattice and extrapolates the results to 
vanishing lattice spacing (a — > 0). A convenient way to carry out this discretisation is to put the fermionic 
variables on the sites of the lattice, whereas the gauge fields are treated as 3 x 3 matrices connecting these sites. 
In this sense, lattice QCD is a classical four-dimensional statistical physics system. One important difference 
compared to three dimensional systems is that the temperature (T) is determined by the additional, Euclidean 
time extension (Nt): T=l/{Nta). Keeping the temperature fixed (such as at the transition point) one can 
reduce a and approach the continuum limit by increasing Nt (see Methods). 

There are several lattice results for the order of the QCD transition (for the two most popular lattice fermion 
formulations see refs and [111), although they have unknown systematics. We emphasise that from the lattice 
point of view two 'ingredients' are necessary to eliminate these systematic uncertainties. 

The first ingredient is to use physical values for the quark masses. Owing to the computational costs this 
is a great challenge in lattice QCD. Previous analyses used computationally less demanding non-physically 
large quark masses. However, these choices have limited relevance. The order of the transition depends on the 
quark mass. For example, in three-fiavour QCD for vanishing quark masses the transition is of first-order. For 
intermediate masses it is most probably a crossover. For infinitely heavy quark masses the transition is again 
first-order. For questions concerning the restoration of chiral symmetry (such as the order of the transition), 
a controlled extrapolation from larger quark masses (such as chiral perturbation theory) is unavailable, and so 
the physical quark masses should be used directly. 

The second ingredient is to remove the uncertainty associated with the lattice discretization. Discretization 
errors disappear in the continuum limit; however, they strongly infiuence the results at non-vanishing lattice 
spacing. In three-flavour unimproved staggered QCD, using a lattice spacing of about 0.28 fm, the flrst-order 
and the crossover regions are separated by a pseudoscalar mass of m^^c^SOO MeV. Studying the same three- 
fiavour theory with the same lattice spacing, but with an improved p4 action (which has different discretization 
errors) we obtain m^_c~70 MeV. In the first approximation, a pseudoscalar mass of 140 MeV (which corresponds 
to the numerical value of the physical pion mass) would be in the first-order transition region, whereas using the 
second approximation, it would be in the crossover region. The different discretisation uncertainties are solely 
responsible for these qualitatively different results [22| . Therefore, the proper approach is to use physical quark 
masses, and to extrapolate to vanishing lattice spacings. Our work eliminates both of the above uncertainties. 

Our goal is to identify the nature of the transition for physical quark masses as we approach the continuum 
limit. We will study the finite size scaling of the lattice chiral susceptibilities ■^t)=5^/(9m^^)(T/F)- log Z, 

where niud is the mass of the light u,d quarks and Ns is the spatial extension. This susceptibility shows 
a pronounced peak around the transition temperature (Tc). For a real phase transition the height of the 
susceptibility peak increases and the width of the peak decreases when we increase the volume. For a first-order 
phase transition the finite size scaling is determined by the geometric dimension, the height is proportional 
to V, and the width is proportional to 1/V. For a second-order transition the singular behaviour is given by 
some power of V, defined by the critical exponents. The picture would be completely different for an analytic 
crossover. There would be no singular behaviour and the susceptibility peak does not get sharper when we 
increase the volume; instead, its height and width will be V independent for large volumes. 

Figure [1] shows the susceptibilities for the light quarks for Nt —A and 6, for which we used aspect ratios 
r = Ns/Nt ranging from 3 to 6 and 3 to 5, respectively. A clear signal for an analytic crossover for both lattice 
spacings can be seen. However, these curves do not say much about the continuum behaviour of the theory. 
In principle a phenomenon as unfortunate as that in the three- flavour theory could occur 22|, in which the 
reduction of the discretization effects changed the nature of the transition for a pseudoscalar mass of «140 MeV. 

Because we are interested in genuine temperature effects we subtract the T=0 susceptibility and study only 
the difference between T^O and T=0 at different lattice spacings. To do it properly, when we approach the 
continuum limit the renormalization of x has to be performed. This leads to rn^Ax, which we study (see 
Methods). 

To give a continuum result for the order of the transition we carry out a flnite size scaling analysis of the 
dimensionless quantity T*/(m^Ax) directly in the continuum limit. For this study we need the height of the 
susceptibility peaks in the continuum limit for flxed physical volumes. The continuum extrapolations are done 
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using four different lattice spacings (A''f=4,6,8 and 10). The volumes at different lattice spacings are fixed in 
units of Tc, and thus VTj?—3'^ ,4'^ and 5^ were chosen. (In three cases the computer architecture did not allow 
us to take the above ratios directly. In these cases, we used the next possible volume and interpolated or 
extrapolated. The height of the peak depends weakly on the volume, so these procedures were always safe.) 
Altogether we used twelve different lattice volumes ranging from 4 • 12"^ to 10 • 48'^ at T > 0. For the T — runs 
lattice volumes from 24 • 12'^ up to 56 • 28"^ were used. The number of trajectories were between 1500 and 8000 
for T > and between 1500 and 3000 for T = 0, respectively. Figure [5] shows the continuum extrapolation for 
the three different physical volumes. The Nt—i results are slightly off but the Nt— 6,8 and 10 results show a 
good a^(xl/Nf scaling. 

Having obtained the continuum values for T"^ / {m? lS.x) at fixed physical volumes, we study the finite size 
scaling of the results. Figure [3] shows our final results. The volume dependence strongly suggests that there is 
no true phase transition but only an analytic crossover in QCD. 

Methods. An introduction to lattice QCD at T—Q and Ty^O can be found in refs [1^ and 24 1 for example. 
The detailed form of our Symanzik improved gauge and stout-link improved staggered fermionic action can 
be found in ref. [l^. (Staggered QCD with Uf^A flavours uses the fourth-root description for the fermionic 



determinant, the relevance of which is recently intensively discussed, see ref. 26| and references therein.) In 
this work we also used a stout-smearing level of 2. Note that stout-link improvement makes the staggered 
fermion taste symmetry violation small even at moderate lattice spacings (for an illustration, see Fig. 1 of 
ref. (25|V Taste symmetry violation errors -as with other discretization errors of the staggered formalism- scale 
as a? and disappear only after extrapolating to the continuum limit. We carried out this extrapolation and 
explicitly checked that taste symmetry violations for our smallest lattice spacings are already within this a? 
scaling regime, and so we conclude that the extrapolation is reliable. 

In previous staggered analyses the gauge configurations were produced by the R-algorithm (27| at a given 
step size. The step size is an intrinsic parameter of the algorithm, which has to be extrapolated to zero. Instead 
of using the approximate R-algorithm, we used [i^l the exact rational hybrid Monte-Carlo algorithm [2^ in 
large scale simulations. This technique, which we apply also in this paper, expresses the fractional powers of 
the Dirac operator by rational functions. 

In our simulations we approach the continuum limit along the line of constant physics (LCP). The LCP is 
defined by relationships between the bare lattice parameters (gauge coupling g and lattice bare quark masses for 
light quarks rUud and strange quark ms). These relationships ensure that the physics (such as ratios of physical 
observables) remain constant, while changing any of the parameters. We note that the LCP is unambiguous 
(independent of the physical quantities, which are used to define the above relationships) only in the vicinity 
of the continuum limit. The present work uses the LCP defined by fixing two ratios to their physical values: 
mi<-/TO7r=3.689 and /i<-/mT=1.185 {rriK is the mass of the kaon, is its decay constant and is the mass 
of the pion). Figure [H shows the LCP used in this work. 

Using the above action, simulation algorithm and parameter set based on our LCP, we performed simulations 
at T^Q and T—Q. At T^Q we always used the physical quark masses at four different temporal extensions: 
iVt=4,6,8 and 10 lattices. The aspect ratios [Ns/Nt) were taken to be 3,4 and 5 (for the Nt — A case also 6), 
which results in a factor-of-8 change in the volume {V). In the chirally broken phase (our zero temperature 
simulations always belong to this class) chiral perturbation theory can be used to extrapolate to the physical 
values of the light quark masses in a controlled manner. Therefore we used four pion masses (rriTr «235, 300, 355 
and 405 MeV), which are somewhat larger than the physical one. The spatial extension was chosen to ensure 
Ns'mTr>4:. The computational requirement for the present work was 0.8 teraflopyears. 

The renormalization of the chiral susceptibility can be done by taking the second derivative of the free 
energy density (/) with respect to the renormalized mass (m^). We apply the usual definition: //T^ = 
-N^[\ogZ{Ns,Nt)/{NtNf) - \ogZ{Nso,Nto)/{NtoN%)]. This quantity has a correct continuum limit. The 
subtraction term is obtained at r==0, for which simulations are carried out on lattices with A^sOi spatial 
and temporal extensions (otherwise at the same parameters of the action). The bare light quark mass is 
related to by the mass renormalisation constant mr=Zm-Tnud- Note that Zm falls out of the combination 
m^9^/9TO^=m^^9^/9m^^. Thus, m^^ [x{Ns,Nt) — x(iVso, iVto)] also has a continuum hmit (for its maximum 
values for different TVt, and in the continuum limit we use the shorthand notation m^Ax). For the subtraction, 
T=0 simulations are performed and chiral perturbation theory is applied. 
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Figure 1: Susceptibilities for the light quarks for Nt=4: (left panel) and for Nt=6 (right panel) as a function 
of where g is the gauge couphng (T grows with The largest volume is eight times bigger than the 

smallest one, so a first-order phase transition would predict a susceptibility peak that is eight times higher (for 
a second-order phase transition the increase would be somewhat less, but still dramatic). Instead of such a 
significant change we do not observe any volume dependence. Error bars are s.e.m. 




Figure 2: Normalised susceptibilities T^/(m^Ax) for the light quarks for aspect ratios r=3 (left panel) r=4 
(middle panel) and r=5 (right panel) as functions of the lattice spacing. Continuum extrapolations are carried 
out for all three physical volumes and the results are given by the leftmost blue diamonds. Error bars are s.e.m 
with systematic estimates. 
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Figure 3: Continuum extrapolated susceptibilities T'^/{m'^Ax) as a function of l/{Tj^V). For true phase tran- 
sitions the infinite volume extrapolation should be consistent with zero, whereas for an analytic crossover the 
infinite volume extrapolation gives a non- vanishing value. The continuum-extrapolated susceptibilities show no 
phase-transition-like volume dependence, though the volume changes by a factor of five. The V^oo extrapo- 
lated value is 22(2) which is llcr away from zero. For illustration, we fit the expected asymptotic behaviour for 
first-order and 0(4) (second order) phase transitions shown by dotted and dashed lines, which results in chance 
probabilities of 10^^^ (7 x 10" ^■^), respectively. Error bars are s.e.m with systematic estimates. 
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Figure 4: The line of constant physics. We show our choice for (strange quark mass) and 20m„d (u,d quark 
masses) in lattice units as functions of 
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